# Create Figures and Tables for Manuscript from Lucid Sample
# Created: 01.03.2019

library(sandwich)
library(lmtest)
library(ri)
library(ggplot2)
library(RColorBrewer)
library(gridExtra)
library(xtable)
library(texreg)
library(dplyr)
library(Hmisc)
library(estimatr)

rm(list = ls())

source(".../02 Lucid Data Clean.R")
source(".../Go-to Stats Functions.R")

figdir <- "PATH TO FOLDER TO OUTPUT FIGURES/TABLES"

######################################
# Functions to Make Figures and Tables
######################################

## Treatment Group Means Plot Function

histmeans <- function(data, treat_var, treat_levels, response, wts, treat_labs, alpha, its, 
                      ylab, title, addc, ymin, ymax) {
  
  tmeans <- boot_stats(data = data, treat_var = treat_var, treat_levels = treat_levels, 
                       response = response, wts = wts, treat_labs = treat_labs, alpha = alpha,
                       its = its)
  
  tmeans$treatlab <- factor(tmeans$treatlab, levels = treat_labs)
  
  p <-  
    ggplot(tmeans, aes(x = treatlab, y = mean, fill = treatlab)) + 
    geom_bar(position = position_dodge(), stat ="identity") +
    geom_errorbar(aes(ymin = cilb, ymax = ciub),
                  width=.2,                    
                  position=position_dodge(.9), colour = "black") +
    
    guides(fill = FALSE) +
    xlab("") + ylab(ylab) + ggtitle(title) + coord_cartesian(ylim = c(ymin, ymax)) +
    
    scale_fill_grey(start = 0.30, end = 0.90) + 
    
    theme(panel.background = element_blank(),
          panel.grid.major = element_blank(), 
          panel.grid.minor = element_blank(),
          axis.text.x = element_text(size = 20, colour = "black"), axis.ticks.x = element_blank(),
          axis.text.y = element_text(size = 20, colour = "black"),
          plot.title = element_text(size = 20, face = "bold", hjust = 0.5), 
          axis.title.y = element_text(size = 20))
  
  for (j in 1:length(treat_levels)) {
    p <- p + annotate("text", x = j, y = tmeans[j,4] + addc, size = 7,
                      label = paste(round(tmeans[j,2], digits = 3),
                                    paste("n = ", tmeans[j,5], sep = ""), 
                                    sep = "\n"))
  }
  
  return(p)
  #return(tmeans)
}

## HMID/UI Treatment Group Differences Table Function

difftable_hmui <- function(data, treat_var, treat_var_c, g, wts) {
  
  require(Hmisc)
  require(dplyr)
  
  if (g == "hm") {
    sampcovs <- c("Support", "Spillover", "Favor Status Quo", "Change Benefit",
                  "Ineq", "Benefits wealthy?", "Knowledge")
    
    resp <- c("_main", "_spill", "_quo", "_s2", 
              "_ineq", "_skew", "_k")
  }
  
  else {
    sampcovs <- c("Support", "Spillover", "Favor Status Quo", "Change Benefit", "Change Premiums",
                  "Ineq", "Benefits wealthy?", "Knowledge")
    
    resp <- c("_main", "_spill", "_quo", "_s2", "_s3", 
              "_ineq", "_skew", "_k")
  }
  
  fill <- rep(NA, length(sampcovs))
  
  sampdif <- data.frame(var = resp, Outcomes = sampcovs, Regressive = fill, Progressive = fill, pval = fill, adj_pval = fill)
  
  attach(data)
  
  for(i in 1:length(resp)) {
    sampdif[sampdif$var == resp[i], "Regressive"] <- wtd.mean(data[treat_var == "reg", paste(g, resp[i], sep = "")], 
                                                              data[treat_var == "reg", wts])
    
    sampdif[sampdif$var == resp[i], "Progressive"] <- wtd.mean(data[treat_var == "pro", paste(g, resp[i], sep = "")], 
                                                               data[treat_var == "pro", wts])
    
    sampdif[sampdif$var == resp[i], "pval"] <- boot_test(data, treat_var_c, c("reg", "pro"), 
                                                         paste(g, resp[i], sep = ""), wts, 0.05, 10000)[3]
  }
  
  detach(data)
  
  # Adjust p-values
  
  if (g == "hm") {
    sampdif[3:4, 6] <- p.adjust(sampdif[3:4, 5], "BH")
    sampdif[5:7, 6] <- p.adjust(sampdif[5:7, 5], "BH")
  }
  
  else {
    sampdif[3:5, 6] <- p.adjust(sampdif[3:5, 5], "BH")
    sampdif[6:8, 6] <- p.adjust(sampdif[6:8, 5], "BH")
  }
  
  sampdif[1:2, 6] <- sampdif[1:2, 5]
  
  sampdif$var <- gsub(pattern = "_", replacement = "", x = sampdif$var)
  colnames(sampdif) <- c("Variable", "Outcome", "Regressive", "Progressive", "p-value", "Adj p-value")
  
  # return(print(xtable(sampdif), include.rownames = FALSE))
  return(sampdif)
}

## Nutrition/Job Training Treatment Group Differences Table Function

difftable_nj <- function(data, treat_var, treat_var_c, g, wts) {
  
  require(Hmisc)
  require(dplyr)
  
  sampcovs <- c("Support", "% Take-up", "Cents/$ Divert Total", "Compare to Alt", "Effect on Taxes",
                "Prog Goal", "Admin", "Govt Misuse", "Citizen Misuse", 
                "less $40k v. $40-75k", "$40-75k v. $75k more",
                "Claim = Need", "Not Claim = No Need", "Effort", "Depend on Govt")
  
  resp <- c("_r", "_tu_4", "_divert", "_alt", "_tax", 
            "_w_1", "_w_2", "_w_3", "_w_4", 
            "_dist_1", "_dist_2", 
            "_need1", "_need2", "_effort", "_depend")
  
  fill <- rep(NA, length(sampcovs))
  
  sampdif <- data.frame(var = resp, Outcomes = sampcovs, Direct = fill, Indirect = fill, pval = fill, adj_pval = fill)
  
  attach(data)
  
  for(i in 1:length(resp)) {
    sampdif[sampdif$var == resp[i], "Direct"] <- wtd.mean(data[treat_var == "direct", paste(g, resp[i], sep = "")], 
                                                          data[treat_var == "direct", wts])
    
    sampdif[sampdif$var == resp[i], "Indirect"] <- wtd.mean(data[treat_var == "indirect", paste(g, resp[i], sep = "")], 
                                                            data[treat_var == "indirect", wts])
    
    sampdif[sampdif$var == resp[i], "pval"] <- boot_test(data, treat_var_c, c("direct", "indirect"), 
                                                         paste(g, resp[i], sep = ""), wts, 0.05, 10000)[3]
  }
  
  detach(data)
  
  # Adjust p-values
  
  sampdif[1, 6] <- sampdif[1, 5]
  sampdif[2:5, 6] <- p.adjust(sampdif[2:5, 5], "BH")
  sampdif[6:9, 6] <- p.adjust(sampdif[6:9, 5], "BH")
  sampdif[10:15, 6] <- p.adjust(sampdif[10:15, 5], "BH")
  
  sampdif$var <- gsub(pattern = "_", replacement = "", x = sampdif$var)
  colnames(sampdif) <- c("Variable", "Outcome", "Direct", "Indirect", "p-value", "Adj p-value")
  
  # return(print(xtable(sampdif), include.rownames = FALSE))
  return(sampdif)
}

## Save Plot Function

outplots <- function(plots, figdir) {
  for (i in 1:length(plots)) {
    pdf(file = paste(figdir, plots[i], ".pdf", sep = ""), width = 10, height = 10)
    print(eval(parse(text = plots[i])))
    dev.off()
  }
}

## Prep Table for Output

preptable <- function(table) {
  require(dplyr)
  table <- select(table, -starts_with("p-value"), -starts_with("Variable"))
  names(table)[4] <- "p-value"
  return(table)
}

########################
# MAIN RESULTS FOR PAPER
########################

# Nutrition and Job Training Support

n_sup_rev <- histmeans(x, "n_treat_v", c("direct", "indirect"), "n_r", 
                        "rwt", c("\nDirect", "\nIndirect"), 0.05, 10000,
                        "[-2,2] = Strongly oppose...Strongly support\n", 
                        "", 
                        0.2, 0, 1.25)

n_sup_rev

j_sup_rev <- histmeans(x, "j_treat_v", c("direct", "indirect"), "j_r", 
                        "rwt", c("\nDirect", "\nIndirect"), 0.05, 10000,
                        "[-2,2] = Strongly oppose...Strongly support\n", 
                        "", 
                        0.2, 0, 1.25)

j_sup_rev

outplots(c("n_sup_rev", "j_sup_rev"), figdir)

boot_test(x, "n_treat_v", c("direct", "indirect"), "n_r", 
          "rwt", 0.05, 10000)

boot_test(x, "j_treat_v", c("direct", "indirect"), "j_r", 
          "rwt", 0.05, 10000)

n_sup_ate <- difference_in_means(n_r ~ n_treat_v, data = x,
                    condition1 = "direct",
                    condition2 = "indirect",
                    weights = rwt)
n_sup_ate

j_sup_ate <- difference_in_means(j_r ~ j_treat_v, data = x,
                    condition1 = "direct",
                    condition2 = "indirect",
                    weights = rwt)

j_sup_ate

n_sup_ate$coefficients / sqrt(wtd.var(x$n_r[x$n_treat_v == "direct"], 
                              weights = x$rwt[x$n_treat_v == "direct"]))

j_sup_ate$coefficients / sqrt(wtd.var(x$j_r[x$j_treat_v == "direct"], 
                              weights = x$rwt[x$j_treat_v == "direct"]))


# Nutrition and Job Training Mechanisms

n_mech_difs <- difftable_nj(x, n_treat_v, "n_treat_v", "n", "rwt")
j_mech_difs <- difftable_nj(x, j_treat_v, "j_treat_v", "j", "rwt")

n_mech_difs

n_mech_rev <- select(n_mech_difs, -Variable, -"p-value")
n_mech_rev      
names(n_mech_rev)[4] <- "p-value"
n_mech_rev

print(xtable(n_mech_rev), only.contents = TRUE, include.rownames = FALSE,
      include.colnames = TRUE, floating = FALSE,
      hline.after = 0,
      file = paste(figdir, "n_mech_rev.tex", sep = ""))

j_mech_difs

j_mech_rev <- select(j_mech_difs, -Variable, -"p-value")
j_mech_rev      
names(j_mech_rev)[4] <- "p-value"
j_mech_rev

print(xtable(j_mech_rev), only.contents = TRUE, include.rownames = FALSE,
      include.colnames = TRUE, floating = FALSE,
      hline.after = 0,
      file = paste(figdir, "j_mech_rev.tex", sep = ""))

n_alt_ate <- difference_in_means(n_alt ~ n_treat_v, data = x,
                                 condition1 = "direct",
                                 condition2 = "indirect",
                                 weights = rwt)
n_alt_ate

j_alt_ate <- difference_in_means(j_alt ~ j_treat_v, data = x,
                                 condition1 = "direct",
                                 condition2 = "indirect",
                                 weights = rwt)
j_alt_ate

##

j_tax_ate <- difference_in_means(j_tax ~ j_treat_v, data = x,
                                 condition1 = "direct",
                                 condition2 = "indirect",
                                 weights = rwt)
j_tax_ate

##

n_dist_1_ate <- difference_in_means(n_dist_1 ~ n_treat_v, data = x,
                                   condition1 = "direct",
                                   condition2 = "indirect",
                                   weights = rwt)
n_dist_1_ate

## 

n_effort_ate <- difference_in_means(n_effort ~ n_treat_v, data = x,
                                    condition1 = "direct",
                                    condition2 = "indirect",
                                    weights = rwt)
n_effort_ate

# HMID and UI Support

hm_sup_rev <- histmeans(x[x$policy == "hm", ], "policy_treat", c("reg", "pro"), "hm_main",
                         "rwt", c("\nRegressive", "\nProgressive"), 0.05, 10000,
                         "[-2,2] = Decreased a lot...Increased a lot\n",
                         "",
                         0.1, 0, 0.8)
hm_sup_rev

ui_sup_rev <- histmeans(x[x$policy == "ui", ], "policy_treat", c("reg", "pro"), "ui_main",
                         "rwt", c("\nRegressive", "\nProgressive"), 0.05, 10000,
                         "[-2,2] = Decreased a lot...Increased a lot\n",
                         "",
                         0.1, 0, 0.8)
ui_sup_rev

outplots(c("hm_sup_rev", "ui_sup_rev"), figdir)

boot_test(x[x$policy == "hm", ], "policy_treat", c("reg", "pro"), "hm_main", 
          "rwt", 0.05, 10000)

boot_test(x[x$policy == "ui", ], "policy_treat", c("reg", "pro"), "ui_main", 
          "rwt", 0.05, 10000)

ui_sup_ate <- difference_in_means(ui_main ~ policy_treat, data = x[x$policy == "ui", ],
                                 condition1 = "reg",
                                 condition2 = "pro",
                                 weights = rwt)

ui_sup_ate

# HMID and UI Mechanisms

hm_mech_difs <- difftable_hmui(x[x$policy == "hm", ], policy_treat, "policy_treat", "hm", "rwt")
ui_mech_difs <- difftable_hmui(x[x$policy == "ui", ], policy_treat, "policy_treat", "ui", "rwt")

hm_mech_difs

hm_mech_rev <- select(hm_mech_difs, -Variable, -"p-value")
hm_mech_rev
names(hm_mech_rev)[4] <- "p-value"
hm_mech_rev

print(xtable(hm_mech_rev), only.contents = TRUE, include.rownames = FALSE,
      include.colnames = TRUE, floating = FALSE,
      hline.after = 0,
      file = paste(figdir, "hm_mech_rev.tex", sep = ""))

ui_mech_difs

ui_mech_rev <- select(ui_mech_difs, -Variable, -"p-value")
ui_mech_rev
names(ui_mech_rev)[4] <- "p-value"
ui_mech_rev

print(xtable(ui_mech_rev), only.contents = TRUE, include.rownames = FALSE,
      include.colnames = TRUE, floating = FALSE,
      hline.after = 0,
      file = paste(figdir, "ui_mech_rev.tex", sep = ""))

#################################
# SUPPLEMENTAL FIGURES AND TABLES
#################################

# Nutrition/Job Training Unweighted Results

n_difs_all_nwts <- difftable_nj(x, n_treat_v, "n_treat_v", "n", NULL)

n_difs_all_nwts

j_difs_all_nwts <- difftable_nj(x, j_treat_v, "j_treat_v", "j", NULL)

j_difs_all_nwts

# HMID/UI Unweighted Results

hm_difs_all_nwts <- difftable_hmui2(x[x$policy == "hm", ], 
                                    policy_treat, "policy_treat", "hm", NULL)
hm_difs_all_nwts

ui_difs_all_nwts <- difftable_hmui2(x[x$policy == "ui", ], 
                                    policy_treat, "policy_treat", "ui", NULL)
ui_difs_all_nwts


n_difs_all_nwts <- preptable(n_difs_all_nwts)
j_difs_all_nwts <- preptable(j_difs_all_nwts)

print(xtable(n_difs_all_nwts), only.contents = TRUE, include.rownames = FALSE,
      include.colnames = TRUE, floating = FALSE,
      hline.after = 0,
      file = paste(figdir, "n_difs_all_nwts", ".tex", sep = ""))

print(xtable(j_difs_all_nwts), only.contents = TRUE, include.rownames = FALSE,
      include.colnames = TRUE, floating = FALSE,
      hline.after = 0,
      file = paste(figdir, "j_difs_all_nwts", ".tex", sep = ""))

hm_difs_all_nwts <- preptable(hm_difs_all_nwts)
ui_difs_all_nwts <- preptable(ui_difs_all_nwts)

print(xtable(hm_difs_all_nwts), only.contents = TRUE, include.rownames = FALSE,
      include.colnames = TRUE, floating = FALSE,
      hline.after = 0,
      file = paste(figdir, "hm_difs_all_nwts", ".tex", sep = ""))

print(xtable(ui_difs_all_nwts), only.contents = TRUE, include.rownames = FALSE,
      include.colnames = TRUE, floating = FALSE,
      hline.after = 0,
      file = paste(figdir, "ui_difs_all_nwts", ".tex", sep = ""))

# Partisan Results

n_difs_dems_ywts <- difftable_nj(x[x$dem == 1, ], n_treat_v, "n_treat_v", "n", "rwt")
n_difs_reps_ywts <- difftable_nj(x[x$rep == 1, ], n_treat_v, "n_treat_v", "n", "rwt")

n_difs_dems_nwts <- difftable_nj(x[x$dem == 1, ], n_treat_v, "n_treat_v", "n", NULL)
n_difs_reps_nwts <- difftable_nj(x[x$rep == 1, ], n_treat_v, "n_treat_v", "n", NULL)

j_difs_dems_ywts <- difftable_nj(x[x$dem == 1, ], j_treat_v, "j_treat_v", "j", "rwt")
j_difs_reps_ywts <- difftable_nj(x[x$rep == 1, ], j_treat_v, "j_treat_v", "j", "rwt")

j_difs_dems_nwts <- difftable_nj(x[x$dem == 1, ], j_treat_v, "j_treat_v", "j", NULL)
j_difs_reps_nwts <- difftable_nj(x[x$rep == 1, ], j_treat_v, "j_treat_v", "j", NULL)


n_difs_dems_ywts <- preptable(n_difs_dems_ywts)
n_difs_reps_ywts <- preptable(n_difs_reps_ywts)
n_difs_dems_nwts <- preptable(n_difs_dems_nwts)
n_difs_reps_nwts <- preptable(n_difs_reps_nwts)

j_difs_dems_ywts <- preptable(j_difs_dems_ywts)
j_difs_reps_ywts <- preptable(j_difs_reps_ywts)
j_difs_dems_nwts <- preptable(j_difs_dems_nwts)
j_difs_reps_nwts <- preptable(j_difs_reps_nwts)

print(xtable(n_difs_dems_ywts), only.contents = TRUE, include.rownames = FALSE,
      include.colnames = TRUE, floating = FALSE,
      hline.after = 0,
      file = paste(figdir, "n_difs_dems_ywts", ".tex", sep = ""))

print(xtable(n_difs_reps_ywts), only.contents = TRUE, include.rownames = FALSE,
      include.colnames = TRUE, floating = FALSE,
      hline.after = 0,
      file = paste(figdir, "n_difs_reps_ywts", ".tex", sep = ""))

print(xtable(n_difs_dems_nwts), only.contents = TRUE, include.rownames = FALSE,
      include.colnames = TRUE, floating = FALSE,
      hline.after = 0,
      file = paste(figdir, "n_difs_dems_nwts", ".tex", sep = ""))

print(xtable(n_difs_reps_nwts), only.contents = TRUE, include.rownames = FALSE,
      include.colnames = TRUE, floating = FALSE,
      hline.after = 0,
      file = paste(figdir, "n_difs_reps_nwts", ".tex", sep = ""))

print(xtable(j_difs_dems_ywts), only.contents = TRUE, include.rownames = FALSE,
      include.colnames = TRUE, floating = FALSE,
      hline.after = 0,
      file = paste(figdir, "j_difs_dems_ywts", ".tex", sep = ""))

print(xtable(j_difs_reps_ywts), only.contents = TRUE, include.rownames = FALSE,
      include.colnames = TRUE, floating = FALSE,
      hline.after = 0,
      file = paste(figdir, "j_difs_reps_ywts", ".tex", sep = ""))

print(xtable(j_difs_dems_nwts), only.contents = TRUE, include.rownames = FALSE,
      include.colnames = TRUE, floating = FALSE,
      hline.after = 0,
      file = paste(figdir, "j_difs_dems_nwts", ".tex", sep = ""))

print(xtable(j_difs_reps_nwts), only.contents = TRUE, include.rownames = FALSE,
      include.colnames = TRUE, floating = FALSE,
      hline.after = 0,
      file = paste(figdir, "j_difs_reps_nwts", ".tex", sep = ""))

###

hm_difs_dems_ywts <- difftable_hmui(x[x$policy == "hm" & x$dem == 1, ], 
                                     policy_treat, "policy_treat", "hm", "rwt")
hm_difs_dems_nwts <- difftable_hmui(x[x$policy == "hm" & x$dem == 1, ], 
                                     policy_treat, "policy_treat", "hm", NULL)

hm_difs_reps_ywts <- difftable_hmui(x[x$policy == "hm" & x$rep == 1, ], 
                                     policy_treat, "policy_treat", "hm", "rwt")
hm_difs_reps_nwts <- difftable_hmui(x[x$policy == "hm" & x$rep == 1, ], 
                                     policy_treat, "policy_treat", "hm", NULL)

ui_difs_dems_ywts <- difftable_hmui(x[x$policy == "ui" & x$dem == 1, ], 
                                     policy_treat, "policy_treat", "ui", "rwt")
ui_difs_dems_nwts <- difftable_hmui(x[x$policy == "ui" & x$dem == 1, ], 
                                     policy_treat, "policy_treat", "ui", NULL)

ui_difs_reps_ywts <- difftable_hmui(x[x$policy == "ui" & x$rep == 1, ], 
                                     policy_treat, "policy_treat", "ui", "rwt")
ui_difs_reps_nwts <- difftable_hmui(x[x$policy == "ui" & x$rep == 1, ], 
                                     policy_treat, "policy_treat", "ui", NULL)

hm_difs_dems_ywts <- preptable(hm_difs_dems_ywts)
hm_difs_dems_nwts <- preptable(hm_difs_dems_nwts)
hm_difs_reps_ywts <- preptable(hm_difs_reps_ywts)
hm_difs_reps_nwts <- preptable(hm_difs_reps_nwts)

ui_difs_dems_ywts <- preptable(ui_difs_dems_ywts)
ui_difs_dems_nwts <- preptable(ui_difs_dems_nwts)
ui_difs_reps_ywts <- preptable(ui_difs_reps_ywts)
ui_difs_reps_nwts <- preptable(ui_difs_reps_nwts)

print(xtable(hm_difs_dems_ywts), only.contents = TRUE, include.rownames = FALSE,
      include.colnames = TRUE, floating = FALSE,
      hline.after = 0,
      file = paste(figdir, "hm_difs_dems_ywts", ".tex", sep = ""))

print(xtable(hm_difs_dems_nwts), only.contents = TRUE, include.rownames = FALSE,
      include.colnames = TRUE, floating = FALSE,
      hline.after = 0,
      file = paste(figdir, "hm_difs_dems_nwts", ".tex", sep = ""))

print(xtable(hm_difs_reps_ywts), only.contents = TRUE, include.rownames = FALSE,
      include.colnames = TRUE, floating = FALSE,
      hline.after = 0,
      file = paste(figdir, "hm_difs_reps_ywts", ".tex", sep = ""))

print(xtable(hm_difs_reps_nwts), only.contents = TRUE, include.rownames = FALSE,
      include.colnames = TRUE, floating = FALSE,
      hline.after = 0,
      file = paste(figdir, "hm_difs_reps_nwts", ".tex", sep = ""))

print(xtable(ui_difs_dems_ywts), only.contents = TRUE, include.rownames = FALSE,
      include.colnames = TRUE, floating = FALSE,
      hline.after = 0,
      file = paste(figdir, "ui_difs_dems_ywts", ".tex", sep = ""))

print(xtable(ui_difs_dems_nwts), only.contents = TRUE, include.rownames = FALSE,
      include.colnames = TRUE, floating = FALSE,
      hline.after = 0,
      file = paste(figdir, "ui_difs_dems_nwts", ".tex", sep = ""))

print(xtable(ui_difs_reps_ywts), only.contents = TRUE, include.rownames = FALSE,
      include.colnames = TRUE, floating = FALSE,
      hline.after = 0,
      file = paste(figdir, "ui_difs_reps_ywts", ".tex", sep = ""))

print(xtable(ui_difs_reps_nwts), only.contents = TRUE, include.rownames = FALSE,
      include.colnames = TRUE, floating = FALSE,
      hline.after = 0,
      file = paste(figdir, "ui_difs_reps_nwts", ".tex", sep = ""))

table(x[x$policy == "hm", "dem"])
table(x[x$policy == "hm", "rep"])

table(x[x$policy == "ui", "dem"])
table(x[x$policy == "ui", "rep"])

# Ideology Results

n_difs_libs_ywts <- difftable_nj(x[x$libs == 1, ], n_treat_v, "n_treat_v", "n", "rwt")
n_difs_cons_ywts <- difftable_nj(x[x$cons == 1, ], n_treat_v, "n_treat_v", "n", "rwt")

n_difs_libs_nwts <- difftable_nj(x[x$libs == 1, ], n_treat_v, "n_treat_v", "n", NULL)
n_difs_cons_nwts <- difftable_nj(x[x$cons == 1, ], n_treat_v, "n_treat_v", "n", NULL)

j_difs_libs_ywts <- difftable_nj(x[x$libs == 1, ], j_treat_v, "j_treat_v", "j", "rwt")
j_difs_cons_ywts <- difftable_nj(x[x$cons == 1, ], j_treat_v, "j_treat_v", "j", "rwt")

j_difs_libs_nwts <- difftable_nj(x[x$libs == 1, ], j_treat_v, "j_treat_v", "j", NULL)
j_difs_cons_nwts <- difftable_nj(x[x$cons == 1, ], j_treat_v, "j_treat_v", "j", NULL)

##

hm_difs_libs_ywts <- difftable_hmui(x[x$policy == "hm" & x$libs == 1, ], 
                                    policy_treat, "policy_treat", "hm", "rwt")
hm_difs_cons_ywts <- difftable_hmui(x[x$policy == "hm" & x$cons == 1, ], 
                                    policy_treat, "policy_treat", "hm", "rwt")

hm_difs_libs_nwts <- difftable_hmui(x[x$policy == "hm" & x$libs == 1, ], 
                                    policy_treat, "policy_treat", "hm", NULL)
hm_difs_cons_nwts <- difftable_hmui(x[x$policy == "hm" & x$cons == 1, ], 
                                    policy_treat, "policy_treat", "hm", NULL)

ui_difs_libs_ywts <- difftable_hmui(x[x$policy == "ui" & x$libs == 1, ], 
                                    policy_treat, "policy_treat", "ui", "rwt")
ui_difs_cons_ywts <- difftable_hmui(x[x$policy == "ui" & x$cons == 1, ], 
                                    policy_treat, "policy_treat", "ui", "rwt")

ui_difs_libs_nwts <- difftable_hmui(x[x$policy == "ui" & x$libs == 1, ], 
                                    policy_treat, "policy_treat", "ui", NULL)
ui_difs_cons_nwts <- difftable_hmui(x[x$policy == "ui" & x$cons == 1, ], 
                                    policy_treat, "policy_treat", "ui", NULL)
##

n_difs_libs_ywts <- preptable(n_difs_libs_ywts)
n_difs_cons_ywts <- preptable(n_difs_cons_ywts)
n_difs_libs_nwts <- preptable(n_difs_libs_nwts)
n_difs_cons_nwts <- preptable(n_difs_cons_nwts)

j_difs_libs_ywts <- preptable(j_difs_libs_ywts)
j_difs_cons_ywts <- preptable(j_difs_cons_ywts)
j_difs_libs_nwts <- preptable(j_difs_libs_nwts)
j_difs_cons_nwts <- preptable(j_difs_cons_nwts)

hm_difs_libs_ywts <- preptable(hm_difs_libs_ywts)
hm_difs_cons_ywts <- preptable(hm_difs_cons_ywts)
hm_difs_libs_nwts <- preptable(hm_difs_libs_nwts)
hm_difs_cons_nwts <- preptable(hm_difs_cons_nwts)

ui_difs_libs_ywts <- preptable(ui_difs_libs_ywts)
ui_difs_cons_ywts <- preptable(ui_difs_cons_ywts)
ui_difs_libs_nwts <- preptable(ui_difs_libs_nwts)
ui_difs_cons_nwts <- preptable(ui_difs_cons_nwts)

table(x$libs)
table(x$cons)

table(x$libs[x$policy == "hm"])
table(x$cons[x$policy == "hm"])

table(x$libs[x$policy == "ui"])
table(x$cons[x$policy == "ui"])

print(xtable(n_difs_libs_ywts), only.contents = TRUE, include.rownames = FALSE,
      include.colnames = TRUE, floating = FALSE,
      hline.after = 0,
      file = paste(figdir, "n_difs_libs_ywts", ".tex", sep = ""))

print(xtable(n_difs_cons_ywts), only.contents = TRUE, include.rownames = FALSE,
      include.colnames = TRUE, floating = FALSE,
      hline.after = 0,
      file = paste(figdir, "n_difs_cons_ywts", ".tex", sep = ""))

print(xtable(n_difs_libs_nwts), only.contents = TRUE, include.rownames = FALSE,
      include.colnames = TRUE, floating = FALSE,
      hline.after = 0,
      file = paste(figdir, "n_difs_libs_nwts", ".tex", sep = ""))

print(xtable(n_difs_cons_nwts), only.contents = TRUE, include.rownames = FALSE,
      include.colnames = TRUE, floating = FALSE,
      hline.after = 0,
      file = paste(figdir, "n_difs_cons_nwts", ".tex", sep = ""))

print(xtable(j_difs_libs_ywts), only.contents = TRUE, include.rownames = FALSE,
      include.colnames = TRUE, floating = FALSE,
      hline.after = 0,
      file = paste(figdir, "j_difs_libs_ywts", ".tex", sep = ""))

print(xtable(j_difs_cons_ywts), only.contents = TRUE, include.rownames = FALSE,
      include.colnames = TRUE, floating = FALSE,
      hline.after = 0,
      file = paste(figdir, "j_difs_cons_ywts", ".tex", sep = ""))

print(xtable(j_difs_libs_nwts), only.contents = TRUE, include.rownames = FALSE,
      include.colnames = TRUE, floating = FALSE,
      hline.after = 0,
      file = paste(figdir, "j_difs_libs_nwts", ".tex", sep = ""))

print(xtable(j_difs_cons_nwts), only.contents = TRUE, include.rownames = FALSE,
      include.colnames = TRUE, floating = FALSE,
      hline.after = 0,
      file = paste(figdir, "j_difs_cons_nwts", ".tex", sep = ""))

print(xtable(hm_difs_libs_ywts), only.contents = TRUE, include.rownames = FALSE,
      include.colnames = TRUE, floating = FALSE,
      hline.after = 0,
      file = paste(figdir, "hm_difs_libs_ywts", ".tex", sep = ""))

print(xtable(hm_difs_cons_ywts), only.contents = TRUE, include.rownames = FALSE,
      include.colnames = TRUE, floating = FALSE,
      hline.after = 0,
      file = paste(figdir, "hm_difs_cons_ywts", ".tex", sep = ""))

print(xtable(hm_difs_libs_nwts), only.contents = TRUE, include.rownames = FALSE,
      include.colnames = TRUE, floating = FALSE,
      hline.after = 0,
      file = paste(figdir, "hm_difs_libs_nwts", ".tex", sep = ""))

print(xtable(hm_difs_cons_nwts), only.contents = TRUE, include.rownames = FALSE,
      include.colnames = TRUE, floating = FALSE,
      hline.after = 0,
      file = paste(figdir, "hm_difs_cons_nwts", ".tex", sep = ""))

print(xtable(ui_difs_libs_ywts), only.contents = TRUE, include.rownames = FALSE,
      include.colnames = TRUE, floating = FALSE,
      hline.after = 0,
      file = paste(figdir, "ui_difs_libs_ywts", ".tex", sep = ""))

print(xtable(ui_difs_cons_ywts), only.contents = TRUE, include.rownames = FALSE,
      include.colnames = TRUE, floating = FALSE,
      hline.after = 0,
      file = paste(figdir, "ui_difs_cons_ywts", ".tex", sep = ""))

print(xtable(ui_difs_libs_nwts), only.contents = TRUE, include.rownames = FALSE,
      include.colnames = TRUE, floating = FALSE,
      hline.after = 0,
      file = paste(figdir, "ui_difs_libs_nwts", ".tex", sep = ""))

print(xtable(ui_difs_cons_nwts), only.contents = TRUE, include.rownames = FALSE,
      include.colnames = TRUE, floating = FALSE,
      hline.after = 0,
      file = paste(figdir, "ui_difs_cons_nwts", ".tex", sep = ""))


# Sample Characteristics

covs <- c("age", "female", "black", "latino", "dem", "rep", "college", "hhimid")
covs_labs <- c("Age", "Female (Yes = 1)", "Black (Yes = 1)", "Latino (Yes = 1)", 
               "Democrats (Yes = 1)", "Republicans (Yes = 1)",
               "College (Yes = 1)", "Income (see notes)")

dl <- data.frame(Covariate = covs_labs, var = covs, 
                 Mean = rep(NA, length(covs)), SD = rep(NA, length(covs)),
                 Mean_WT = rep(NA, length(covs)), SD_WT = rep(NA, length(covs)))

for (i in 1:length(covs)) {
  dl[i,3] <- mean(x[ ,covs[i]], na.rm = T)
  dl[i,5] <- wtd.mean(x[ ,covs[i]], x$rwt, na.rm = T)
  
  if (dl[i,2] %in% c("age", "hhimid")) {
    dl[i,4] <- sd(x[ ,covs[i]], na.rm = T)
    dl[i,6] <- sqrt(wtd.var(x[ ,covs[i]], x$rwt, na.rm = T))
  }
}

print(xtable(select(dl, -var)), only.contents = TRUE, include.rownames = FALSE,
      include.colnames = TRUE, floating = FALSE,
      hline.after = 0,
      file = paste(figdir, "lucid_demos", ".tex", sep = ""))
